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ABSTRACT 

At redshifts z <; 2, most of the baryons reside in the smooth intergalactic medium 
which is responsible for the low column density Lya forest. This photoheated gas fol- 
lows a tight temperature-density relation which introduces a cut-off in the distribution 
of widths of the Lya absorption lines (6-parameters) as a function of column density. 
We have measured this cut-off in a sample of nine high resolution, high signal-to-noise 
quasar spectra, and determined the thermal evolution of the intergalactic medium in 
the redshift range 2.0-4.5. At redshift z ~ 3, the temperature at the mean density 
shows a peak and the gas becomes nearly isothermal. We interpret this as evidence 
for the reionization of He II. 

Key words: cosmology: miscellaneous - intergalactic medium - quasars: absorption 
lines 



1 INTRODUCTION 

According to the standard big bang model, the primordial 
hydrogtpii and helium cumpiisiiig the hile-igalacLie medium 



1996). Measurements of the Hell Lya op acity suggest that 
helium may have rei o nized around z ~ 3 ( Heap_e£_aL_200C ; 



(IGM) was hot and highly ionized at early times. As the 

universe expanded, the hot plasma cooled adiabatically, be- 
coming almost completely neutral at a redshift of z ~ 10 3 . 
The IGM remained neutral until the first stars and quasars 
began to produce ionizing photons. Eventually, the ioniz- 
ing radiation became intense enough to reionize hydrogen 
and later, because of its higher ionization potential, to fully 
reionize helium. Since the thermal evolution of the IGM de- 
pends strongly on its reionization history, it can be used as 
a probe of the end of the 'dark ages' of cosmic history, when 



Refiners et al. 1997 ; lakobsen et al. 1994 ; Davidsen et ah 
1996j ; [Anderson et al. 1999| ). This would fit in with evidence 



for a hardening of the UV background around this time, as 
derived from the ratio of Siiv/Civ in h i gh redshift qu asar 
absorption lines ( Songaila fc Cowie 1996] ; |5ongaila 1998) ) , al- 
though both the obser vational result and its interpretation 
are still controversial ( Boksenberg, Sargent & Rauch 199£ ; 
Giroux fc Shull 1997|). 



Rees 1994; 


Huifc Gnedin 1997; 


Haehnelt & Steinmetz 1998|). 


The at 


>sence of Gunn- Peterson absorption ( Gunn & Pe- 



terson 1965) in quasar spectra, i.e. the complete absorption 
of quas ar light blueward of the H I and He H Lva wavelengths 
require ! that hydrogen must have been highly ionized by 



The resonant Lya absorption by residual low levels of 
neutral hydrogen along the line of sight to a quasar pro- 
duces a forest of absorption lines. Although many of the 
basic observational facts about the Lya forest at high red- 
shift (z ~ 2-5) had been established before the 10 m tele- 
scope era, the advent of the Keck telescope has lead to 
much larger data samples at much higher signal-to-noise ra- 



tio than hitherto available (e.g. Hu et al. 199£ 



199£) and helium by 



5 (Schneider, Schmidt & Gunn 1991; pongaila et al 



Lu et al 



199(i ; Kirkman fc Tytler 1997[ ). The observational progress 



2.5 (Davidsen, Kriss & Zheng 



has bee n matched on the theoret ical side by semi-analytic 
models (e.g. Bi fc Davidsen 1997) and cosmological hydro- 



man 1991 



The lata presented herein were obtained at the W. M. Keck 
Observatory, which is operated as a scientific partnership among 
the California Institute of Technology, the University of Califor- 
nia and the National Aeronautics and Space Administration. The 
Observatory was made possible by the generous financial support 
of the W. M. Keck Foundation. 



simulation s (s.g. Cen et al. 199 4 ; Zhang. Anninos & Nor 



lYtitjcau, Mucket fc Kates 1995^ Hernquist et al 



1996; Miralda-Escude et al. 1996) which together with the 



new data are now beginning to yield significant quantita- 
tive cosmological constraints (see Rauch 1998 for a recent 
review) . 

These simulations show that the low column density 
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(N Jj 10 14 ' 5 cm -2 ) absorption lines arise in a smoothly vary- 
ing IGM of low density contrast (S ^ 10), which contains 
most of the baryons in the universe. Since the overden- 
sity is only mildly non-linear, the physical processes gov- 
erning this medium are well understood and relatively easy 
to model. On large scales the dynamics are determined by 
gravity, while on small scales gas pressure is important. 
Since shock heating is unimportant for the low-density gas, 
the interplay between photoionization heating and adia- 
batic cooling due to the expansion of the universe results 
in a tight temperature-density relation, which is well de- 
scribed by a powe r-law for densities ar ound the cosmic mean, 
T = T (p/p) 7_1 ( |Hui fc Gnedin 1997| ). This relation is gen- 
erally referred to as the 'equation of state' (even though the 
true equation of state is that of an ideal gas). 

For models with abrupt reionization, the IGM becomes 
nearly isothermal (7 ~ 1) at the redshift of reionization. Af- 
ter reionization, the temperature at the mean density (To) 
decreases while the slope (7 — 1) increases because higher 
density regions undergo increased photoheating and expand 
less rapidly. Eventually, the imprints of the reionization his- 
tory are washed out and the equation of state approaches 
an asymptotic state, 7 = 1.62, Tp oc [Q.th 2 / ^H^ h 2 ^ 1 ' 



(Miralda-Escude & Rees 1994; Hui & Gnedin 1997; rheum 



et al. 1998). However, the timescale for recombination cool- 
ing in the low density IGM is never small compared to the 
age of the universe for z ^ 20 and inverse Compton cooling 
of free electrons off the cosmic microwave background is only 
efficient for z <; 5. Consequently, unless both hydrogen and 
helium were fully reionized at redshifts considerably higher 
than this, the gas will have retained some memory of when 
and how it was reionized. 

A standard way of analyzing Lya forest spectra is to 
decompose them into a set of distinct absorption lines, as- 
sumed to have Voigt profiles (e.g. Carswell et al. 1987). Vari- 
ous broadening mechanisms, such as Hubble broadening (the 
differential Hubble flow across the absorber), peculiar and 
thermal velocities contribute to the line widths (Meiksin 
1994; Hui & Rutledge 1999; Theuns, Schaye & Haehnelt 
2000). However, there exists a lower limit to the line- width, 
set by the temperature of the gas. Because the physical den- 
sity of the IGM correlates strongly with the column density 
of the absorption lines, this results in a cut-off in the distri- 
bution of line widths (b-parameters) as a function of column 
density, which traces the equation of state of the gas (Schaye 
et al. 1999, hereafter STLE; Ricotti, Gnedin & Shull 2000; 
Bryan & Machacek 2000). Hence we can infer the equation 
of state of the IGM by measuring the minimum Lya line 
width as a function of column density. 

Here, we measure the b(N) cut-off in nine high resolu- 
tion, high S/N quasar spectra, spanning the redshift range 
2.0-4.5. We use hydrodynamic simulations to calibrate the 
relations between the parameters of the b(N) cut-off and the 
equation of state. By applying these relations to the obser- 
vations, we are able to measure the evolution of the equation 
of state over the observed redshift range. We find that the 
thermal evolution of the IGM is drastically different from 
that predicted by current models. The temperature peaks 
at 2 ~ 3, which, together with supporting evidence from 
measurements of the Hell opacity and the Siiv/Civ ratios, 
we interpret as evidence for the second reionization of he- 



Table 1. Quasar spectra used 

QSO -em 

Q1100-264 2.14 

Q2343+123 2.52 

Q1442+293 2.67 

Q1107+485 3.00 

Q1425+604 3.20 

Q1422+231 3.62 

APM 08279+5255 3.91 

Q0000-262 4.11 

Q2237-061 4.55 



Hum (Hen — > Hem). Ricotti et al. ( 2000 ) recently applied 
a similar technique to published lists of Voigt profile fits. A 
comparison with the method and results of Ricotti et al. is 
given in section ^. 

This paper is organized as follows. In sections |^ and 
|§] we describe the observations and the simulations respec- 
tively. We discuss the difference between evolution of the 
fe-distribution and evolution of the temperature in section ^j. 
In section [| we briefly describe our method for measuring 
the equation of state, before we present our results in sec- 
tion ^| Systematic errors are discussed in section g Finally, 
we discuss and summarize the main results in section H. 



2 OBSERVATIONS 

We analyzed a sample of nine quasar spectra, spanning the 
redshift range z cul = 2.14-4.55 (Table |l|). The spectra of 
Q1100-264 and APM 08279+5255 were kindly provided by 
R. Carswell and S. Ellison respectively. All spectra were 
taken with the high-resolution spectrograph (HIRES, Vogt 
et al. 1994) on the Keck telescope, except the spectrum of 
Q1100— 264, which was taken with the UCL echelle spec- 
trograph of the Anglo Australian Telescope. Details on the 
data and reduction procedures, as we ll as t he continuum fit- 
ting, can be fo und in Carswell et al. fll99l| ) for Q1100-264, 
Ellison et al. (|l999| ) for APM 0827 9+525 5 and Barlow & 
Sargent (1997) and Rauch et al. (1997) for the others. 
The nominal velocity resolution (FWHM) was 8kms _1 for 
Q1100— 264 and 6.6kms _1 for the others and the data were 
rebinned onto 0.04 A pixels on a linear wavelength scale. The 
signal to noise ratio per pixel is typically about 50, except 
for Q1100-264 for which it is about 20. 



In order to avoid confusion with the Ly/3 forest, only 
the region of a spectrum between the quasar's Ly/3 and Lya 
emission lines was considered. In addition, spectral regions 
close to the quasar (typically 8-10 Mpc, but 21 h" 1 Mpc 
for APM 08279+5255 and 32 h" 1 Mpc for Q1100-264) were 
omitted to avoid proximity effects. Regions thought to be 
contaminated by metals and damped Lya lines were re- 
moved (metal line regions were identified by correlating with 
metal lines redwards of the quasar's Ly-alpha emission line 
and with strong Hi lines). The absorption features in the 
remaining spectral regions were fitted with V oigt profiles 
using the same auto mated version of VPFIT (Webb 1987; 
Carswell et al. 1987) as was used for the simulated spectra. 
Using a fully automatic fitting program invariably results 
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in a few 'bad fits'. However, given that there is no unique 
way of decomposing intrinsically non-Voigt absorption lines 
into a set of discrete Voigt profiles, it is essential to apply 
the same algorithm to simulated and observed spectra. Since 
'bad fits' will also occur in the synthetic spectra, we have 
made no attempt to correct them. 

The Lya forest of a single quasar spans a considerable 
redshift range (Az ~ 0.5). In order to minimize the effects 
of redshift evolution and S /N variation across a single spec- 
trum, we divided each Lya forest spectrum into two parts 
of equal length. STLE showed that their algorithm for mea- 
suring the cut-off of the b(N) distribution is relatively in- 
sensitive to the number of absorption lines, the statistical 
variance is almost the same for e.g. 150 and 300 lines. Hence 
little information is lost by analyzing narrow redshift bins 
if the absorption line density is high. The two halves of the 
spectra were analyzed separately and each was compared 
with its own set of simulated spectra (see section ^|). For the 
two lowest redshift quasars the number of absorption lines 
is too small to split the data in half. Hence each quasar, 
except for Q1100— 264 and Q2343+123, provides two nearly 
independent data sets. 

The complete set of absorption line samples is listed in 
Table ^. The median, minimum and maximum redshifts of 
the absorption lines used to determine the b(N) cut-off are 
listed in columns 2-4. The minimum column density con- 
sidered was set to 10 12 ' 5 cm -2 for all samples, since blends 
dominate at lower column densities. The maximum column 
densities are listed in column 5, they were determined by 
the following considerations: (a) the cut-off can be measured 
more accurately if the column density interval is larger; (b) 
we only want to measure the cut-off for column densities 
that correspond to the density range for which the gas fol- 
lows a power-law temperature-density relation. The total 
number of absorption lines in each sample is listed in col- 
umn 6 (only lines for which VPFIT gives relative errors in 
both the b-parameter and column density less than 0.25 are 
considered) . Finally, column 6 lists the pivot column density 
of the power-law cut-off, b = &at (iV/7Vo) r_1 , that was fit to 
the data (c.f. section [5]). 



Table 2. Observed absorption line samples. Columns 2—4 cor- 
respond to the median, minimum and maximum redshifts of 
the absorption lines in the sample. Column 5 is the maximum 
column density considered (the minimum column density is al- 
ways 10 12 ' 5 cm -2 ) and column 6 contains the total number of 
absorption lines used to determine the b(N) cut-off. The last 
column is the pivot column density of the power-law cut-off, 
b = b No {N/N ) T - 1 . 



Sample 


z 




^max 


log TVmax 


# lines 


log iVo 


1100 


1.96 


1.85 


2.09 


14.0 


44 


13.0 


2343 


2.29 


2.23 


2.38 


14.2 


48 


13.0 


1442a 


2.21 


2.10 


2.33 


14.5 


87 


13.0 


1442b 


2.50 


2.33 


2.63 


14.5 


90 


13.0 


1107a 


2.59 


2.48 


2.71 


14.5 


76 


13.0 


1107b 


2.84 


2.71 


2.95 


14.5 


91 


13.0 


1425a 


2.66 


2.55 


2.88 


14.5 


94 


13.0 


1425b 


3.00 


2.89 


3.14 


14.5 


118 


13.0 


1422a 


3.08 


2.91 


3.22 


14.5 


145 


13.0 


1422b 


3.37 


3.22 


3.53 


14.5 


152 


13.0 


0827a 


3.23 


3.15 


3.42 


14.5 


105 


13.4 


0827b 


3.55 


3.43 


3.70 


14.5 


106 


13.4 


0000a 


3.72 


3.61 


3.81 


14.8 


77 


13.5 


0000b 


3.91 


3.81 


4.01 


14.8 


97 


13.5 


2237a 


3.84 


3.69 


4.02 


14.8 


140 


13.5 


2237b 


4.31 


4.15 


4.43 


14.8 


127 


13.5 




4.5 


4.0 


3.5 


z 

3.0 


2.5 


2.0 



0.0 
-0.2 
-0.4 
-0.6 
-0.8 
-1.0 



O 2237 

• 0000 
□ 0827 
■ 1422 
O 1425 

♦ 1107 
A 1442 
A 2343 
x 1100 



-1.2 1 ...... 

0.75 0.70 



0.65 0.60 0.55 
log(l+z) 



0.50 



0.45 



Scatter plots of the 6(iV)-distribution for all observed 
samples (Table ^) are shown in Fig. |l|. The solid lines are 
the measured cut-offs, the vertical dashed lines indicate the 
maximum column density used for fitting the cut-off. There 
are clear differences between the samples. We will show in 
section ^| that even a non-evolving ^-distribution would im- 
ply a strong thermal evolution. Several samples contain a 
few lines that fall far below the cut-off. These lines, which 
have no significant effect on the measured cut-off, are most 
likely blends or unidentified metal lines. 

Fig. ^ shows the effective optical depth, r e fi = 
— \n({F)), of the observed samples as a function of (decreas- 
ing) redshift. The ionizing background in the simulations 
was rescaled to match these effective optical depths. The 
scatter is small, considering that most data points repre- 
sent iust half of a Lya forest spectrum. Note that Rauch et 
al. ( 1997 ) studied the opacity of the forest using seven out of 
nine of the quasars from this sample. They found a slightly 
less rapid increase with redshift, because they rebinned the 
data into three redshift bins, centered on z = 2, 3 and 4. 



Figure 2. The observed effective optical depth as a function of 
redshift. The dashed line, which has a slope of 4.5, is drawn to 
guide the eye. Horizontal error bars indicate redshift intervals 
used, vertical error bars are 1 a errors as determined from a boot- 
strap analysis using chunks of 4 A. 



3 SIMULATIONS 

In order to calibrate the relation between the parameters of 
the b(N) cut-off and the effective equation of state, we have 
simulated eight variants of the currently favoured flat, scale 
invariant, cosmological constant dominated cold dark mat- 
ter model, which vary only in their heating rates (Table ^J). 
The calibration was repeated for each of the observed sam- 
ples of absorption lines listed in Table |^. Synthetic spectra 
were computed along 1200 random lines of sight through the 
simulation box at the nearest redshift output (Az = 0.25). 
The background flux was rescaled such that the mean effec- 
tive optical depth in the simulated spectra matches that of 
the observed sample. Each spectrum was convolved with a 
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Figure 1. b(7V)-Distributions for the observed samples listed in Table ^ Crosses indicate positions of absorption lines, errors are not 
displayed. Only lines that are used for the determination of the cut-off are shown, i.e. lines for which VPFIT gives relative errors in both 
b and N smaller than 25 per cent. Solid lines are the measured cut-offs. Vertical dashed lines indicate the maximum column density used 
when fitting the cut-off. Horizontal dashed lines arc identical and correspond to the cut-off of sample 1425b. Note that evolution of the 
b(N) cut-off cannot be interpreted as evolution of the equation of state in any straightforward way, because the density-column density 
relation changes with redshift (see section^). Furthermore, the cut-offs in different panels can only be compared directly if they contain 
a similar number of absorption lines. 



Gaussian with full width at half maximum (FWHM) identi- 
cal to that of the observations and resampled onto pixels of 
the same size. The noise properties of the observed spectrum 
were computed as a function of flux and imposed on the sim- 
ulated spectra. The resulting s pectr a were continuum fitted 
as described in Theuns et al. (1998). Finally, Voigt profiles 
were fitted using the same automated version of VPFIT as 
was used for the observed spectra. We will refer to the sam- 
ple of lines drawn from the synthetic spectra of simulation 
X, designed to mimic the observed spectrum Y as as model 
X-Y, e.g. model Ll-1442a. 

All models have a total matter density Q m = 0.3, vac- 
uum energy density Qa = 0.7, baryon density Qth 2 = 0.019, 
present day Hubble constant Ho = 65kms _1 Mpc _1 and 
the amplitude of the initial power spectrum is normalized 
to as — 0.9. The IGM is assumed to be of primordial com- 
position with a helium abundance of 0.24 by mass and is 
photoionized and photoheated by the UV-background from 
quasars. The exact cosmology and UV-background is unim- 



Table 3. Simulations used for calibrating the relation between 
the b(N) cut-off and the effective equation of state. 



Model 




Hx 


Comment 


L0.3 


1/3 





low To 


LI 


1 





reference model 


L2 


2 





high T 


L3 


3 





very high T 


Lx 


1 


1 


low 7, high To 


Lx2.5 


1 


5/2 


very low 7, high To 


Lx5 


1 


5 


very low 7, very high To 


Lie 


1 





high 7 



portant for this analysis, as is the normalization (see sec- 
tion Q). Although these parameters may affect the equation 
of state of the IGM in the simulations, they do not change 
its relation to the b(N) cut-off. 
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The numerical simulations used in this paper follow the 
evolution of a periodic, cubic region of the univ erse and are 



perfor med with a modifie d version of HYDRA (Couchman 
Thoma > fc Pearce 1995|) , which uses smooth partic le hy- 
drodynamics ( [Lucy 1970] ; |Gingold fc Monaghan 1977| ). The 



simulations employ 64 J gas particles and 64 a cold dark mat- 
ter particles in a box of of comoving size 3.85 Mpc, so the 
particle mass is 1.14 x 10 6 M for the gas and 6.51 x 10 6 M 
for the dark matter. 

Our reference model, LI, is photoionized and photo- 
heated by the UV-background from quasars as computed 
by Haardt & Madau (1996, hereafter HM), using the opti- 
cally thin limit. Models L0.3, L2 and L3 are identical, except 
that we have multiplied the helium photoheating rates (col- 
umn tHo in Table fa) by factors of 1/3, 2 and 3 respectively 
(keeping the ionization rates constant). The effective helium 
photoheating rate may be higher than computed in the opti- 
cally thin limit because of radiative transfer effects (e.g. Abel 
& Haehnelt 1999). Model Lx is identical to model LI, ex- 
cept that we have included Compton heating by the h ard X - 
ray background as computed by Madau & Efstathiou ( 199!: ) 
(column in Table [|) . For a highly ionized plasma, the en- 
ergy input per particle from Compton scattering of free elec- 
trons is independent of the density. Hence, Compton heating 
tends to flatten the effective equation of state. We have used 
this fact to artificially construct models with low values of 
7, by multiplying the X-ray heating rates by (unrealistic) 
factors of 2.5 and 5 for models Lx2.5 and Lx5 respectively. 
Finally, model Lie is identical to model LI, except that we 
have set the ionization and heating rates for H I and He I for 
redshifts between 6 and 10 equal to those at z = 6. In this 
model Hi and Hel ionize early (at z = 10), which drives 7 
to larger values. 

In addition to the models listed in Table ^, we have per- 
formed some simulations to investigate possible systematic 
effects. We simulated model LI twice with lower resolutions 
(2 x 54 3 and 2 x 44 3 particles) and model L3 in a larger box, 
but with the same resolution (5 ft -1 Mpc and 2 x 128 3 parti- 
cles instead of 2.5 h~ x Mpc and 2 x 64 3 particles). Model LI 
was simulated twice more with lower normalizations of the 
initial power spectrum (cr% = 0.65 and 0.4 instead of 0.9). 



EVOLUTION OF THE 6-DISTRIBUTION VS 
THERMAL EVOLUTION 



Williger et al. ( |l994j ) and Lu et al. ( |l996| ) found that the 
^-parameters at z ~ 4 are smaller than at z = 2-3. Kim 
et al. (1997) showed that the increase in the line widths 
with decreasing redshift continues over the range z = 3.5 to 
2.1. It is tempting to interpret these results as evidence for 
an increase in the temperature To with decreasing redshift. 
However, we will show in this section that the 6- values are 
smaller at higher redshift even for models in which To is 
higher, as is the case for models in which the universe is 
fully reionized by z = 4. 

As pointed out by STLE, any statistic that is sensitive 
to the temperature of the absorbing gas, will in general de- 
pend on both the amplitude, To, and the slope, 7, of the 
equation of state. This is because temperature is a function 
of density and the absorbing gas is, in general, not all at 
the mean density of the universe. After reionization, To will 
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Figure 3. The temperature at a given density contrast as a 
function of redshift for model LI, which uses the HM ionizing 
background. Although the temperature at the mean density, To, 
increases with redshift (solid line) , the temperature at slight over- 
densities is almost constant, or even decreasing with redshift. 



decrease and 7 will increase with time. Consequently, the 
evolution of the temperature at a given overdensity can be 
very different from the evolution of To. This is illustrated 
in Fig. where the temperature Ts at a density contrast 
S = p/p — lis plotted as a function of redshift for model LI. 
Even though the temperature at the mean density decreases 
with time (solid line), the temperature at a density contrast 
as little as 2 remains almost constant. 

The general expansion of the universe ensures that the 
column density corresponding to a fixed overdensity is a 
strongly increasing function of redshift. In fact, most of the 
evolution of the Lya forest can be understood in terms of 
the resulting scaling of the optical depth (e.g. Hernquist et 
al. 1996; Dave et al. 1999; Machacek et al. 1999). When in- 
terpreting the evolution of the ^-distribution, one therefore 
has to keep in mind that: (a) at fixed column density, ab- 
sorption lines at higher redshift will correspond to absorbers 
of smaller overdensities; (b) the evolution of the tempera- 
ture at a fixed overdensity depends on the evolution of both 
To and 7. Together these effects can conspire to make the 
6-parameters smaller at higher redshift, even when the tem- 
perature To is higher. Fig. rA shows that this will happen for 
models in which the IGM is fully reionized at the observed 
redshifts (z <,4). 

The discussion in this section shows that to derive the 
evolution of To using a statistic that is sensitive to the tem- 
perature of the absorbing gas, one needs to determine the 
evolution of: (1) the temperature of the gas; (2) the over- 
density of the gas and (3) the slope of the equation of state. 
We will see that the uncertainty in 7 is the limiting factor. 



5 METHOD 

STLE demonstrated that the observed cut-off in the dis- 
tribution of 6-parameters as a function of column density 
can be used to measure the equation of state of the IGM. 
In particular, they showed that the b(N) cut-off can be fit- 
ted by a power-law, b = fcjv (iV/iVo) r-1 , whose parameters 
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Figure 4. The b(N) cut-off at a fixed column density 
(10 13 ' 6 cm -2 ) as a function of redshift for model LI, which uses 
the HM ionizing background. The error bars enclose 68 per cent 
confidence intervals around the medians, as determined from 1000 
sets of 300 absorption lines. The cut-off shifts to larger b- values 
with decreasing redshift, even though the temperature at the 
mean density, To, decreases (solid line in Fig. H). 



log &jv and F — 1 are proportional to log T S ( Nq ) and 7 — 1 re- 
spectively. For each observed sample, we first calibrate these 
relations using the simulations and then use them to convert 
the observed cut-offs into measurements of the equation of 
state. 

For each observed sample of absorption lines (Table ^|) 
we go through the following procedure. First mock spectra 
are generated from the 8 simulations listed in Table pi The 
synthetic spectra are processed to give them the same char- 
acteristics (mean absorption, resolution, pixel size and noise 
properties) as the corresponding observed spectra. These are 
then fitted with Voigt profiles using the same automated fit- 
ting package that was used for the observations. For each of 
the eight simulated sets of absorption lines, the b(N) cut-off 
is fitted using the iterative procedure developed by STLE. 
We then use these 8 simulations to calibrate the relations be- 
tween the parameters of the b(N) cut-off, (6jv ,r) and the 
parameters of the equation of state, (T S (n ) , 7)- We mea- 
sure the density contrast corresponding to the pivot col- 
umn density, No, by using the fact logT^(jv ) oc log6jv 
(STLE), essentially because both thermal broadening and 
Jeans smoothing scale as the square root of the tempera- 
ture. 

Fig. ^ illustrates our method for measuring Tg and 
8 (No). In the left panel the intercept of the cut-off, mea- 
sured at No = 10 1 cm -2 in this example, is plotted as 
a function of log To for the simulations of sample 1422a 
(filled circles). As expected, the relation between the inter- 
cept, log 614.0, and log To is not very tight. The large scatter 
arises because log 614.0 is not proportional to log To, but to 
logTa(jv ), where 6 (No) is the density contrast correspond- 
ing to the pivot column density iVo = 10 14 '°cm -2 . Indeed, 
if we plot log 614.0 as a function of logTs = i.e (right panel), 
the scatter becomes very small, implying that 6 (No) ~ 1.6. 

Changing the value of 5 shifts the data points in the 
plot horizontally, the amount depending on the value of 7 in 
the simulation corresponding to the data point. Because we 



do not know a priori what the density contrast correspond- 
ing to No is, we vary 8 and see for which value the scatter 
is minimal. The inset in the right panel shows the total \ 2 
of the linear least squares fit to the data points as a func- 
tion of the density contrast. For this example the scatter is 
minimal for 8 = 1.6, which is the density contrast used in 
the right panel. Using the optically depth weighted density 
- column density relation, introduced by STLE, we find that 
the column density TV = 10 140 cm 2 does indeed correspond 
to a density contrast of about 1.6. 

The dot-dashed lines in Fig. |E| show the 6- value expected 
for pure thermal broadening, 6 = yj2ksT jm v , where ks is 
the Boltzmann constant and m p is the proton mass. The 
relation between log 614.0 and logTj = i.6 lies slightly above 
the pure thermal broadening line, but has the same slope. 
This implies that the widths at the cut-off are dominated by 
thermal broadening, with an additional component that also 
scales as T 1 / 2 . We identify this last component as the dif- 
ferential Hubble flow across the absorber, whose size is set 
by the Jeans smoothing scale and does indeed depend on 
the square root of the temperature. We consistently found 
that minimizing the scatter in the log 6jv -log T$ relation by 
varying S results in a relation that has a slope of about 
0.5 and that the value of 8 found agrees well with direct 
measurements of the density contrast corresponding to the 
column density No- We therefore use this procedure to esti- 
mate 8 and logT^ and conservatively estimate the error in 
the density contrast to be cr(log(l + 8)) — 0.15. Although for 
this work the error from the determination of the density - 
column density relation does not contribute significantly to 
the total error in To, it may well become the limiting factor 
when a larger sample of quasars is used. 

Having measured Tg, 8 and 7 we can compute To, 

logT = logT a -( 7 -l)log(l + 5). (1) 

Each measurement of log6jv and V comes with associated 
errors, which are determined from the bootstrap distribution 
(c.f. STLE). These errors can be converted directly into er- 
rors in logTa and 7 respectively, using the linear relations 
between the cut-off and the equation of state, determined 
from the simulations. To these errors we add in quadrature 
the residual scatter of the data points around the linear fit. 
The error in log To is then given by, 

A 2 (logT ) = A 2 (logT,) + [A(7)log(l+<5)] 2 + 

[( 7 -l)A(log(l + 5))] 2 . (2) 

After having measured the thermal evolution of the 
IGM, we ran a simulation designed to match the observa- 
tions (dashed lines in Fig. ^). The open square in Fig. ^ 
corresponds to this simulation. The difference between the 
evolution in the calibrating simulations and the true evolu- 
tion is particularly large at z ~ 3 (c.f. Fig. |^a), the redshift 
of model 1422a. The fact that the open square follows the 
same relation as the other data points, confirms that a model 
whose equation of state matches the one determined from 
the observations using the methods described in this section, 
does indeed have the observed b(N) cut-off. 
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Figure 5. The intercept of the b(N) cut-off as a function of the temperature at the mean density (left panel) and at a density contrast of 
8 = 1.6 (right panel). Filled circles are the simulations corresponding to the observed sample 1422a. The open square is for a simulation 
whose thermal evolution matches the observations (dashed line in Fig. ^). The dashed line is the least-squares fit for the filled circles. 
The dot-dashed line indicates the fe-value corresponding to pure thermal broadening: b = (2fcsT/m p ) 1//2 . The inset in the right panel 
shows how the (total) \ 2 of the fit varies as a function of 8; \ 2 ls minimum for 8 = 1.6, which is the density contrast corresponding 
to the column density of the intercept of the b(N) cut-off. (10 14 cm - 2 in this example). The difference between the dashed and the 
dot-dashed lines in the right panel is due to the contribution of baryon smoothing to the widths of the lines around the b(N) cut-off. 
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Figure 6. The temperature at the mean density (a) and the the slope of the effective equation of state (b) as a function of redshift. 
Horizontal error bars indicate the redshift interval spanned by the absorption lines, vertical error bars are 1 a errors. The solid lines are 
for our reference model, LI, which uses the HM ionizing background. The dashed lines are for a simulation that was designed to fit the 
data. In this model, which has a much smaller contribution from quasars at high redshift, Hen reionizes at z ~ 3.2. Radiative transfer 
effects on the temperature of the IGM were modeled schematically by increasing the photoheating rates for an optically thin gas during 
reionization by a factor 4 for H I and He II and a factor 2 for He I. Although it is clear that the temperature peaks at z ~ 3 and that the 
gas becomes close to isothermal (7 ~ 1), the present constraints are not sufficient to distinguish between a sharp rise (as indicated by 
the dashed line) and a more gradual increase. 



6 RESULTS 

The measured evolution of the temperature at the mean 
density and the slope of the effective equation of state are 
plotted in Fig. ^| From z ~ 4 to z ~ 3, To increases and the 
gas becomes close to isothermal (7 ~ 1.0). This behavior 
differs drastically from that predicted by models in which 
helium is fully reionized at higher redshift. For example, the 
solid curves correspond to our reference model, LI, which 
uses a uniform metagalactic UV-background from quasars 



as computed by HM and which assumes the gas to be opti- 
cally thin. In this simulation, both hydrogen and helium are 
fully reionized by z ~ 4.5 and the temperature of the IGM 
declines slowly as the universe expands. Such a model can 
clearly not account for the peak in the temperature at 2 ~ 3 
(reduced x 2 f° r the solid curves are 6.8 for To and 3.6 for 
7). Instead, we associate the peak in To and the low value 
of 7 with reheating due to the second reionization of helium 
(Hen — > He 111). 

If reionization of Hell happens locally on a timescale 
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that is short compared to the recombination timescale, 
which for Hem is of the order of the age of the universe 
at z ~ 3, then the energy density injected by photoioniza- 
tion will be proportional to the gas density. Consequently, 
the temperature increase will be independent of the den- 
sity and the equation of state of the IGM will become more 
isothermal. The change in the slope of the equation of state 
at z ~ 3 is thus physically consistent with our interpretation 
of the peak in the temperature at the same redshift. 

The dashed lines in Fig. ^ are for a model that was con- 
structed to fit the data (reduced x 2 is 0.22 for To and 1.38 
for 7). This model, for which stellar sources ionize Hi and 
He 1 by z ~ 5 and quasars ionize He n at z ~ 3.2, has a much 
softer spectrum at high redshift. Before reionization, when 
the gas is optically thick to ionizing photons, the mean en- 
ergy per ph otoionization is much h igher than in the optically 
thin limit (Abel & Haehnelt 199£). We have approximated 
this effect in this simulation by enhancing the photoheating 
rates during reionization, so raising the temperature of the 
IGM. 

Since the simulation assumes a uniform ionizing back- 
ground, the temperature has to increase abruptly (i.e. much 
faster than the gas can recombine) in order to make 7 as 
small as observed. In reality, the low-density gas may be 
reionized by harder photons, which will be the first ionizing 
photons to escape from the dense regions surrounding the 
sources. This would lead to a larger temperature increase in 
the more dilute, cooler regions, resulting in a decrease of 7 
even for a more gradual reionization. Furthermore, although 
reionization may proceed fast locally (as in our small simula- 
tion box), it may be patchy and take some time to complete. 
Hence the steep temperature jump indicated by the dashed 
line, although compatible with the data, should be regarded 
as illustrative only. The globally averaged To could well in- 
crease more gradually which would also be consistent with 
the data. 

Note that if the reionization is patchy, it would give 
rise to large spatial fluctuations in the temperature. Be- 
cause we measure the temperature-density relation from the 
lower cut-off of the 6(iV)-distribution, our results should 
be regarded as lower limits to the average temperature. 
Absorption lines arising in a local, hot ionization bubble 
would not necessarily raise the observed cut-off in the b(N)- 
distribution. 

The errors in Fig. |^ can be directly traced back to the 
corresponding b(/V)-distributions (Fig. uh . Take for example 
0827a, which has an extremely low value of 7, with a very 
large error. This is clearly due to the gap in the 6-distribution 
at log TV ~ 12.5-13.0. The lack of lines in that region could 
be a statistical fluctuation, or it could be an indication of 
large variations in the temperature of the IGM. 



7 SYSTEM ATICS 

In this section we will investigate whether there are any 
systematic effects that could affect our results. 



7.1 Numerical resolution 

The fe-distribution h as been shown to be very sensitive to nu 
merical resolution ( Theuns et al. 1998 ; Bryan et al. 1999 ) 



It is therefore important to check that the lower limit to the 
line widths in our simulations is not set by the numerical res- 
olution. We have resimulated our reference model, LI, which 
is the second coldest model, twice more at a lower resolution. 
The resolution was decreased by decreasing the number of 
particles from 2 x 64 3 to 2 x 54 3 and 2 x 44 3 respectively, 
while keeping the size of the simulation box constant. The 
resulting probability distributions for the intercept and the 
slope of the b(N) cut-off are plotted in Fig. [j] for our lowest 
and highest redshift samples. Only for the lowest resolu- 
tion simulation do the differences become noticeable. The 
intercept increases slightly and the slope becomes slightly 
shallower, indicating that the lines at the low column den- 
sity end are not resolved. We conclude that the simulations 
used for this work have sufficient resolution to provide an 
accurate determination of the effective equation of state of 
the IGM. 



7.2 Simulation box size 

In order to investigate the effect of the simulation box size, 
we resimulated model L3 using a larger box, but with the 
same resolution (5 h~ x Mpc and 2 x 128 3 particles instead of 
2.5 IT X Mpc and 2 x 64 3 particles). The effect of increasing 
the size of the simulation box is more difficult to determine 
than the effect of numerical resolution. When the box size is 
increased, the gas becomes slightly hotter (To increases by 
a few per cent), presumably because shock heating is more 
effective due to the larger infall velocities. This effect is small 
and since it does not affect the relation between the b(N) 
cut-off and the equation of state, it is unimportant for this 
work. We find that using the larger simulation box results 
in higher values of To. The effect is negligible at z ~ 4 and 
increases to AlogT » 0.05 (AT /T « 0.12) at z ~ 2. 
The change in the derived value of 7 is very small (^ 0.05). 
Although To and 7 may change a bit more if we increase the 
box size further, the small difference between the two box 
sizes indicates that the effect of the box size is insignificant 
compared to the statistical errors. 



7.3 Cosmology 

STLE showed that the relation between the cut-off in the 
b(N) distribution and the equation of state is independent of 
the assumed cosmology. However, the initial power spectra 
of the models investigated by STLE were all normalized to 
match the observed abund ance of galaxy clusters at z = 0. 
Bryan & Machacek (2000) claimed that the b-distribution 



depends strongly on the amplitude of the po wer sp ectrum, as 
predicted by t he mo del of Hui & Rutledge (199£). However, 
Theuns et al. (2000) found the dependence to be very weak, 
provided that the line fitting is done using an algorithm, 
like VPFIT, that attempts to deblend absorption lines into 
a set of thermally broadened components. Although the ab- 
sorption features do become broader for models with less 
small-scale power, the curvature in the line centers does not 
change much. Consequently, the total number of Voigt pro- 
file components used by VPFIT will generally increase, but 
the fits to the line centers will change very little. 

To quantify the effect of decreasing the amount of small- 
scale power, we resimulated our reference model LI twice 
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Figure 7. The effect of numerical resolution on the b(N) cut-off. Shown are the probability distributions for the intercept (left) and 
slope (right) of the cut-off in the b(7V)-distributions of models L1-Q1100 (z = 1.97) (top) and Ll-Q2237b (z = 4.27) (bottom), using 
2 X 64 3 (solid), 2 X 54 3 (dashed) and 2 X 44 3 (dot-dashed) particles respectively. The small difference between the intermediate and 
high resolution simulations indicates that the latter has converged. The intercept was measured at a column density of 10 14 cm -2 for 
Ll-Q2237b and 10 13 4 cm" 2 for L1-Q1100. 



using a lower normalization of the initial power spectrum. 
We find that normalizing to as = 0.65 instead of as = 0.9, 
changes the derived values of To by less than 3 per cent for all 
redshifts. The change in 7 is never greater than 0.1. For the 
extreme case of erg = 0.4, the derived values of To are about 
15 per cent lower and 7 differs by about 0.15. We conclude 
that the effect of the uncertainty in the normalization of the 
primordial power spectrum is small. 

7.4 Continuum fitting and the mean absorption 

Errors in the continuum fit of the observed spectra, will 
lead to errors in the effective optical depth. Underestimating 
the observed continuum will decrease the measured effective 
optical depth. Decreasing the mean absorption in the sim- 
ulations will increase the density corresponding to a given 
column density. Hence the slope of the cut-off will remain un- 
changed, but the intercept will increase, although the effect 
is small (STLE). Increasing the intercepts of the calibrat- 
ing simulations will decrease the derived temperature Tg. 
The higher derived density contrast will work in the same 
direction (provided that 7 > 1), resulting in a lower To. 

The measured effective optical depth for sample 2343 
seems to be relatively high compared to the other samples 
(Fig. |^) . We tried scaling the synthetic spectra to the effec- 
tive optical depth corresponding to the dashed line in Fig. ^ 
at the redshift of 2343, which is about 30 per cent lower than 
the measured value. This resulted in an increase of log To by 
0.04 (AT0/T0 w 0.09), while leaving 7 unchanged. 

In addition to errors in the continuum fit of the observed 
spectra, the continuum fitting of the synthetic spectra could 
also lead to systematic errors. We checked this by repeating 
the analysis of the simulations corresponding to our highest 



redshift sample, 2237b, but this time without continuum 
fitting the synthetic spectra. In this case the derived value 
of To would be 9 per cent lower, while the value of 7 would 
be higher by 0.07. Hence systematic effects in the continuum 
fitting of the observed and the synthetic spectra are unlikely 
to be important. 



8 SUMMARY AND DISCUSSION 

We have measured the cut-off in the distribution of line 
widths (b) as a function of column density (N) in a set 
of nine high-quality Lycv forest spectra, spanning the red- 
shift range 2.0-4.5. We emphasized that the evolution of 
the temperature of the intergalactic medium (IGM) cannot 
be derived directly from the evolution of the b(N) distribu- 
tion, the decrease of the overdensity corresponding to a fixed 
column density with redshift has to be taken into account. 
We therefore used hydrodynamic simulations to calibrate 
the relations between the b(N) cut-off and the temperature- 
density relation of the low density gas. The calibration was 
done separately for each observed spectrum, using synthetic 
spectra that were processed to give them identical charac- 
teristics as the observed spectrum. Crucially, Voigt profiles 
were fitted to the real and simulated spectra using the same 
automated fitting package (a modified version of VPFIT). 
We have checked possible systematic errors arising from the 
finite numerical resolution, the finite size of the simulation 
box, the amplitude of the initial power spectrum and contin- 
uum fitting errors in both synthetic and observed spectra. In 
all cases, the effects were small compared to the statistical 
errors. 

The measured thermal evolution differs drastically from 
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the scenario predicted by current models of the ionizing 
background from quasars, in which helium is fully reion- 
ized by z ~ 4.5. The temperature at the mean density, To, 
increases from z ~ 4 to z ~ 3, after which it decreases again 
(Fig. ^). The slope of the equation of state reaches a mini- 
mum at z ~ 3, where it becomes close to isothermal. More 
data at 2 J; 3 is needed to determine whether the rise in 
To is sharp or gradual. These results suggest that the low 
density IGM was reheated from z ~ 4-3, which we interpret 
as reheating associated with the reionization of He II. 

These results are in qualitati ve ag reement with those re- 
ported recently by Ricotti et al. (2000). They used a method 
which relies on the assumption that only thermal broadening 
contributes to the line widths of the absorption lines at the 
peak of the 6-distribution and used approximate simulation 
techniq ues to determine the density - column density rela- 



Ll, which uses the HM ionizing background (solid lines in 
Fig. |6ij), with the data clearly shows that the model under- 
estimates the temperature at z ~ 3. Since the above men- 
tioned mechanisms for increasing the temperature do not 
change the overall shape of the thermal evolution, a change 
in the reionization history is required to bring the data and 
simulations into agreement. Furthermore, the photoheating 
rates around reionization need to be enhanced to account for 
the fact that heating by photons with energies significantly 
above the ionizatio n potential is important when the gas is 
not optically thin ( Abel fc Haehnelt 1999 ), as is generally 
assumed in cosmological simulations. 

There are two other lines of evidence for late reioniza- 
tion of He II. The first are direct measurements of the optical 
depth from He n Lya abs orption, which have s o far been ob- 
tained fo r four quasars (Jakobsen et al. 199^; Davidsen et 



tion. applying their method to published lists of Voigt al. 199E ; Reimers et al. 1997; Anderson et al. 1999; Heap 



profile : its they found that at z ~ 3, 7 is smaller than would 



be expected if the reionization of helium had been completed 
at high redshift. Ricotti et al. have only three data points in 
the range z = 2-4 which all overlap at the 0.5 sigma level, so 
we cannot compare the shape of the temperature evolution. 
However, it is interesting that they measure a temperature 
at 2 = 2 and z = 4 that is about seventy per cent higher 
than is reported here, although their error bars are suffi- 
ciently large to agree with our results at the la level. Since 
Jeans smoothing contributes to the line widths, as Fig. ^ 
indicates (see also Theuns et al. 2000), their method should 
lead to an overestimate of the temperature (by about sixty 
per cen t for the example in Fig. ^[). The contribution of Jeans 



et al. 2000). These observations already provide strong ev- 
idence for a drop in the mean absorption from z ~ 3.0 to 
2.5. The second piece of evidence concerns a change in the 
spectral shape of the ionizing background. As Hell is ion- 
ized, the mean free path of hard UV photons will increase 
and the spectrum o f the UV background will become harder. 
Songaila & Cowie (199(3) and Songaila (1998) have reported 



a rapid increase with decreasing redshift of the Silv/ClV 
ratio at z ~ 3, which they interpreted as evidence for asud- 
den reionization of Hell. However, Boksenberg et al. ( 1996 ) 
found only a gradual change with redshift. The interpre- 
tation of this metal line ratio is complicat ed because local 



stella r radiation is likely to be important (Giroux & Shull 



smooth ng to the line width may be larger for the lines near 
the peak of the ^-distribution than indicated in Fig. [5], which 
is for the narrower lines near the cut-off. 

Although the reionization of helium may not be the 
only process which can explain the peak in To at z ~ 3, 
it appears to be the only process that can simultaneously 
account for the observed decrease in 7. Galactic winds for 
example, would be less important in the low density regions, 
and would therefore result in an increase in the slope of the 
equation of state. Another possible explanation is a hard- 
ening of the ionizing background at z ~ 3, as might be 
expect ed because of the increase in the number of quasars. 



1997p . 

It should be kept in mind that these three different types 
of observations probe different physical structures. Our re- 
sults apply to density fluctuations around the cosmic mean, 
the effective optical depth depends mostly on the neutral 
fraction in the voids and the metal line ratios probe the high 
density peaks. These structures will probably not be ionized 
simultaneously. After the ionization front breaks through 
the haloes surrounding the source of He 11 ionizing photons, 
e.g. a quasar, it will propagate quickly into the voids. The 
filaments, where the recombination rate is much hi gher, will 
get ionized more slowly, starting from the outsid e ( |Miralda~ 



If heliu n had already been ionized, this would still raise the Escude, Haehnelt & Rees 2000, Gnedin 2000). The IGM 



temperature somewhat, but the effect would be stronger in 
the high density regions where the gas recombines faster. 
Hence this would also lead to an increase in 7, contrary to 
what is observed. 



Rec ently, it was shown (Theuns et al. 1998; Bryan et al i 
199S |) t| iat high resolution simulations of the standard cold 
dark model, using the ionizing background computed by HM 
produce a larger fraction of narrow lines than observed at 
z ~ 3. Different authors have prop osed different solutions 
to this problem. Theuns et al. (1998) suggested that the gas 
tem perat ure in the simulations was too low, while Bryan et 
al. (1999) argued that the amplitude of primordial fluctua- 



tions was too high. For a given reionization history, the tem- 
perature in the simulations could for example be increased 
by inc reasing the baryon density and the age of the uni- 
verse ( Theuns et al. 199S ), including Compton heatin g by 
the hard X-ray background (Madau fc Efstathiou 199E) and 



possibl y photo-elect ric heating by dust grains (Nath, Seth 



Shch ekinov 1999 ). A comparison of our reference model 



will still be thick to Lya photons when only a small neutral 
fraction remains in the voids. Hence a drop in the He 11 Lya 
optical depth at z ~ 3 would suggest that the reionization 
of the voids, which co ver most of the volume, is complete 
( Miralda-Escude 1996 ) . It is important to note that the peak 
in To at 2 ~ 3 does not imply that the temperature of the 
general IGM reaches a maximum. In fact, our results imply 
that the temperature of slightly overdense gas (5 <; 2) is 
almost constant because the slope of the equation of state, 
7, is minimum when To is maximum. 

Detailed modeling, probably in the form of large hydro- 
dynamical simulations, which include radiative transfer, is 
required to see whether the various observational constraints 
can be fit into a consistent picture. However, the ingredients 
necessary to explain our discovery of a peak in the tempera- 
ture of the low density IGM at 2 ~ 3 are clear even from our 
crude model: a softer background at high redshift to delay 
helium reionization and enhanced heating rates compared to 
the optically thin limit. Once the evolution of helium heat- 
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ing is understood, the measurements of the temperature at 
higher redshifts can be used to constrain the epoch of hy- 
drogen reionization. 

Finally, we would like to note that because of their 
hard spectrum, quasars tend to ionize helium shortly af- 
ter hydrogen, although the delay depends on the dumpiness 
of the IGM (iMadau, Haardt & Rees 19981). It may there- 



fore be difficult to postpone the reionization of helium until 
z ~ 3, if quasars were responsible for reionizing hydrogen at 
2 > 5. Hence the mounting evidence for helium reionization 
at z ~ 3 suggests that hydrogen was reionized by stars. 
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